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The abundance of (not just) dark matter haloes 
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ABSTRACT 

We study the effect of baryons on the abundance of structures and substructures in 
a ACDM cosmology using a pair of high resolution cosmological simulations from 
the GIMIC project. Both simulations use identical initial conditions, but while one 
contains only dark matter, the other also includes baryons. We find that gas pressure, 
reionisation, supernova feedback, stripping, and truncated accretion systematically 
reduce the total mass and the abundance of structures below ~ 10 12 Mq compared to 
the pure dark matter simulation. Taking this into account and adopting an appropriate 
detection threshold lowers the abundance of observed galaxies with maximum circular 
velocities v max < 100 kms -1 , significantly reducing the reported discrepancy between 
ACDM and the measured HI velocity function of the ALFALFA survey. We also show 
that the stcllar-to-total mass ratios of galaxies with stellar masses of ~ 10 5 — 1O 7 M0 
inferred from abundance matching of the (sub)halo mass function to the observed 
galaxy mass function increase by a factor of ~ 2. In addition, we find that an important 
fraction of low-mass subhaloes are completely devoid of stars. Accounting for the 
presence of dark subhaloes below 1O 1O M0 further reduces the abundance of observable 
objects, and leads to an additional increase in the inferred stellar-to-total mass ratio 
by factors of 2-10 for galaxies in haloes of 10 9 — 10 10 Mq. This largely reconciles 
the abundance matching results with the kinematics of individual dwarf galaxies in 
ACDM. We propose approximate corrections to the masses of objects derived from 
pure dark matter calculations to account for baryonic effects. 

Key words: cosmology: theory - galaxies: formation - galaxies: evolution - galaxies: 
luminosity function, mass function - methods: N-body simulations 



1 INTRODUCTION 

The large-scale evolution of the Universe is determined by 
gravity which drives hierarchical structure formation in an 
expanding space. While the nature of dark matter remains 
an open question at a fundamental level, gravity makes 
no distinction between baryonic and dark matter. Conse- 
quently, numerical simulations of cosmic structure formation 
are generally purely gravitational, or "Dark Matter Only" 
(DMO), which greatly reduces the computational cost and 
complexity compared to gas-dynamical simulations. Over 
the past decades, DMO simulations have played a key part 
in establishing the currently preferr ed ACDM paradigm (e.g. 
iDavis et all 19851 : iFrenk et alll988h and the hierarchical for- 
mation scenario (e.g. IWhite fc Reesl 1 19781 : IWhite fc Frenkl 
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Il99ll ). explaining the observed abundance and correlation 
of structures at different redshifts and over many mass and 
length scales. 

On dwarf galaxy scales (haloes up to 10 11 Mq or u max < 
100 kms -1 , or galaxies with stellar masses up to ~ 1O 9 M0), 
however, cold dark matter simulations appear to overpre- 
dict the amount of structure, to a degree that appears in- 
compatible with observations: the abundance of haloes with 
a maximum circular velocity u max ~ 35 kms -1 exceeds the 
abundance mea sured in the ALFALFA HI survey by a fac - 
tor of ~ 10 (e.g. lZavala et al.ll2009l : iPapastergis et a.l.ll201lh : 
the simulations fail to match the dynamics of satellites, (e.g. 
Parry et al]|2009l ; IStrigari et al]|2010l : iBovlan-Kolchin et all 



Wang 2012); and the stellar-to-total 
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mass ratios inferred from the abunda nces do not match 
those of individua l dwarf galaxies (e.g. ISawala et alj I20T1I : 
iFerrero etalll2012i ). In this work we show that the inclusion 
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of baryons and their relevant astrophysical processes changes 
the abundance of (sub) haloes sufficiently to reconcile CDM 
with observations. 

While the observable matter in the Universe traces the 
underlying total mass distributi on, galaxies are "biased" 
with respect to the dark matter (|Kaiseri[l98i : [5 avis et al.l 
GUI), and the stellar-to-total mass ratio varies greatly 
with galaxy mass, evolutionary stage, and environment. 
Linking simulated structures to observations of galaxies 
must therefore be done indirectly, via methods such as 
(sub) ha lo abundance ma tc hing ( e.g. Vale fc Ostrikerl 



2009: IConrov fc Wechsleil 120091 : iMoster et al.l l201p| : 



Guo et al] |2010|). halo occupation d i stribution models 



(e.g iBenson et all l200d; ISeliakl |2000| ; IScoccimarro et al.1 
l200ll : iBerlind fc Weinberd l2002h . or semi-analytica l mod- 
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Kauffmann et al.l [f 993: 
Cole et al. 19941; ISomerville fc Primackllf 999 ; ICroton et al.l 
20061 : I Bower et alJ 120061 ). While these approaches vary 
greatly in complexity, aim and predictive power, they 
share the fundamental assumption that the total mass 
distribution in the Universe is governed by gravity alone, 
and can be represented by dark matter. 

In this work, we investigate the impact of baryon 
physics on structure formation on smaller scales. We find 
that baryons systematically affect the masses and abun- 
dances of objects in three ways: 

• Interstellar gas is expelled from haloes in simulations 
that include supernova feedback, while intergalactic gas is 
prevented from accreting in the presence of ionising radia- 
tion. 

• Ram-pressure stripping preferentially removes gas from 
low-mass satellite galaxies, and satellites are more easily dis- 
rupted. 

• The shallower potential well caused by the loss of 
baryons subsequently leads to diminished accretion, both 
of baryons and of dark matter. 

Of these, the first two effects can substantially reduce 
the mass of baryons in small haloes, and lower the total 
mass of an object by up to the universal baryon fraction, 
compared to the corresponding DMO simulation. The third 
effect can reduce an object's mass even further, as haloes 
with a lower mass and shallower potential accrete less mate- 
rial. All baryon effects show a strong mass dependence, and 
stripping also depends on environment. 

To quantify the influence of baryons, we use a set of 
simulations of structure formation in identical cosmological 
volumes: a pure dark matter simulation, and one which in- 
cludes a gas physics model that incorporates hydrodynam- 
ics, radiative heating and cooling, a cosmic UV background, 
star formation, chemical evolution and supernova feedback. 
This simulation, part of the "Galaxie s-Intergalactic Me dium 
Interaction Calculation" ( "GIMIC" . ICrain et al"1l2009l ). has 
sufficient dynamic range and volume to resolve objects with 
total mass between 10 9 — 1O 14 M0, from dwarf galaxies to 
massive groups and clusters. It has been used previously 
to study many aspects of galaxy formation , such as en- 
viro nmental depende ncies (|Crain et al.l 120091 ') , disk galax- 
ies ([Sales et al.ll2012l ) and their X-ray coronae (|Crain et al.l 
|2010|). the origin of the baryonic Tully Fisher relation 



Here, our focus is not on the formation of galaxies, but 
on the formation of their dark matter hosts, and in partic- 
ular, on the abundance of haloes and subhaloes of different 
masses and circular velocities, in different environments, and 
over time. Matching individual objects across both simula- 
tions also enables us to understand the mechanisms by which 
baryons affect the abundance of structures, and to derive an 
approximate expression that allows us to correct for bary- 
onic effects in generic DMO simulations. 

The effect of baryons on the structure of haloes has 
been studied extensively, both using analytical mod els (e.g. 
iBlumenthal et al.lll986l ; ISellwood fc McGaugbl l 2005l) and in 
numerical simulations (e .g. lAbadi et al.l |2010| ; iDuffv et al.l 
l20ld : iBrvan et all |2012| ). Whereas dissipa tion universally 
increases the conce ntration of haloes (e.g. iMo et all l20ld : 
iGnedin et al.ll201ll ). a non-adiabatic response to rapid out- 
flows induced by supernovae can er ase the centra l cusp 
in the profiles of low-mass haloes ([Navarro et al.l 1 19961 ; 
iGovernato et al.ll2010l ; iPontzen fc Governatoll2012l ). poten- 
tially resolving the discrepancy between cusped profiles 
predicted by CDM, a nd the cored profiles indic ated by 
some observations ( e.g. | Walker fc Penarrubia|[201ll . but see 
IStrigari et aUl2010l ; IWclf fc BuUockll2012l ). 

Comparisons of the masses and abundances of haloes 
in DMO and gas simulations have a lso been pe r forme d 
previously, by Weinberg et al l (120081); IRudd et all (|200Sl ); 
iDolag et all (|2009l ) and ICui et al.l ~(|2012[ ). but these have 
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focuse d on scales of 10 12 Mq and above. ISemboloni et all 
(|201ll ) and Ivan Daalen et ail (|201ll ) have investigated the 
effects of baryons on estimates of cosmological parameters, 
for k = 10/i/Mpc and larger. On these scales, (|Dolag et al.l 
12009? ) concluded that baryons only have a minor effect on 
halo masses, with adiabatic co ntraction balanc ing feedback, 
while IRudd et al.1 (|2008f ) and (|Cui et al.l [20121) both found 
an increase in the mass of haloes, with ICui et al.l (|2012T ) 
reporting a net increase of 1 — 2% in M200 for haloes of 
10 13,5 M Q and above, and IRudd et al.1 (|2008l ) report an in- 
crease in the cumulative halo mass function of ~ 10% above 
1O 12 /i _1 M0. An opposite effect o f a reduction in the mat - 
ter power spectrum was found by Ivan Daalen et al.l (120111), 
in sim ulations which include AGN feedback. Simha et al.l 
(|2012l ) have compared hydrodynamic simulations to results 
of (sub)halo abundance matching, and found the results to 
be consistent. 

By comparison to these works, our study is based on 
much higher resolution simulations, which allows us to inves- 
tigate the effect of baryons in lower mass haloes. In showing 
that the effect of baryons is highly mass-dependent, we find 
that while our results are consistent with previous results 
for objects of 10 12 Mq and abov^H on smaller scales, baryons 
significantly reduce the mass of haloes and subhaloes, with 
consequences for abundance matching. 

This paper is organised as follows: in Section [2j we de- 
scribe the simulations, and outline the procedures employed 
to identify and link substructures across the two simulations. 
Section [3] contains our main results: in Section |3. II we com- 
pare individual structures in both simulations, and we give 
an analytic expression for the change in the median subhalo 
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Figure 1. Evolution of a slice of 16fe _1 Mpc in depth through the total mass distribution within the high-resolution regions of the DMO 
simulation (top) and GIMIC simulation (bottom), at z = 6, 1, 0.5, and (left to right). Shown in each panel is a scale bar of constant 
comoving length 10h~ 1 Mpc, and its corresponding length in physical coordinates. On large scales, the mass distributions in the DMO 
and GIMIC simulations appear identical. 



mass in Section [32] In Section 13.31 we describe the effect of 
baryons on the mass functions of haloes and subhaloes, while 
in Section f3. 41 we examine the t> ma x function, and compare 
our results to observational HI data. In Section [331 we show 
the effect of baryons on the satellite mass functions of groups 
and clusters. Section 13.61 focuses on dark subhaloes, which 
do not contain any gas or stars. In Section [31 we describe 
how the baryonic effects change the results of abundance 
matching, and compare the inferred stellar-to-total mass ra- 
tios to observations. We summarise our results in Section [5] 
Parameters for the analytical mass correction are given in 
Appendix [A] and numerical convergence is discussed in Ap- 
pendix [B] 



The simulations were performed using the TreePM-SPH 
code P-Gadge t-3, an extens ion of the publicly available 
code Gadget-2 (|Springej|2005h . The cosmological parame- 
ters are identical to those of the Millennium simulations; 
density parameters Q m = 0.25, J1a = 0.75, Qb = 0.045, 
Hubble parameter h — 0.73, power spectrum normalisation 
erg = 0.9, and spectral index n = 1. The evolution of the 
dark matter distribution in the full, high-resolution region 
in both simulations is shown in Fig. [1] By construction, the 
simulated volume is approximately spherical at z = 1.5, and 
becomes more irregular at both higher and lower redshifts. 
On scales of megaparsecs and above, baryons do not signifi- 
cantly affect structure formation, and no differences can be 
seen in the density distribution shown in Fig. [T] 



2 METHODS 

2.1 The Simulations 

Our comparison is based on a set of N-Body simulations 
of the same cosmological volume: a hydrodynamical simu- 
lation ("GIMIC"), per formed by the Virg o Consortium and 
described in detail by ICrain et al.l (|2009h . and a matching 
dark matter only simulation ( "DMO" ) where the total mass 
is composed of dark matter. The simul ations refine a sph er- 
ical Lagrangian volume (labelled Oct in lCrain et al.ll2~009l ) of 
radius ~ 18/i _1 Mpc (a volume of ~ 63 Gpc 3 at z = 0) 
and median density , selec ted from the Millennium Simula- 
tion (|Springel et al.ll2005h . Within the zoom region, addi- 
tional small s cale modes a r e add ed using the zoom method 
described by iPower et all ([20031 1. while the external large 
scale density field is represented by lower resolution dark 
matter particles. 



2.1.1 Initial Conditions 

Both simulations evolve the same number, 2.75 x 10 s , of dark 
matter particles, from z = 127 to z = 0, and are analysed 
at 58 co-temporal snapshots. In the GIMIC simulation, the 
high-resolution region additionally contains an equal num- 
ber of gas particles that can be converted to star particles. 

The initial particle masses are m g = 1.98 x 1O 6 M0 and 
TtinM = 9.05 x 1O 6 M0 for gas and dark matter, respec- 
tively in the GIMIC simulation, and ttidm = 1-98 + 9.05 = 
11.03 x 10 6 M Q in the DMO simulation. In the GIMIC sim- 
ulation, gas particles can be converted into star particles of 
equal mass, and mass can be exchanged between gas and 
star particles in the form of stellar winds and supernovae. 
The gravitational softening scale for all particle types in the 
high resolution region is initially fixed in comoving coordi- 
nates, and chosen so that it reaches a value of 0.5h~ kpc 
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Figure 2. Comparison of halo mass estimators M200,mean (blue), 
Mtophat (red), and M2oo,crit (black), as a function of Mp p for 
all haloes at z = 0. Solid lines show the results of the GIMIC 
simulation, dashed lines show results of the DMO simulation. As 
expected, for all objects, M2oo,mcan gives the highest mass per 
halo, followed by Mj op h a t and M2oo,crit- Importantly, there is no 
systematic bias in mass estimators between the GIMIC or the 
DMO simulation. 




Figure 3. Top panel: cumulative fraction of subhaloes in the 
DMO simulation that are uniquely matched to subhaloes in the 
GIMIC simulation at z = (solid) and z = 1 (dashed), for cen- 
trals (red), satellites (blue), and all subhaloes (black). Bottom 
panel: fraction of all subhaloes at z = that are central in the 
DMO simulation (black), and in the GIMIC simulation (red). In 
general, the fraction of subhaloes that arc matched increases with 
mass, and it is slightly lower at z = compared to jar = 1, par- 
ticularly for satellites. The lower matching fraction among satel- 
lites and its decrease over time is attributable to a higher rate 
of mergers, and a higher probability for divergent evolution be- 
tween the two simulations. Overall, more than 90% of subhaloes 
above 1O 9 M are matched at z = 0. In both simulations, centrals 
account for ~ 80% of subhaloes in the mass range 1O 9 -1O 12 M0. 



Plummer equivalent in physical coordinates at z = 3, where- 
upon it remains fixed in physical coordinates, contracting 
in comoving coordinates as 1/a. To test convergence, both 
simulations were also repeated at lower resolution, with all 
particle masses increased by a factor of eight, and the soften- 
ing increased by a factor of two. Unless we specifically refer 
to the lower resolution simulations for the purpose of com- 
parison, all results are based on the high-resolution simula- 
tions. We discuss convergence in Appendix [B] and consider 
the abundance of (sub)haloes to be unaffected by resolution 
down to 10 9 M Q . 

2.1.2 Physical Model 

While the DMO simulation includes only gravitational in- 
teractions, the GIMIC simulation also includes astrophysical 
processes, such as gas dynamics, photo-ionisation, radiative 
element-by-element cooling, star formation and supernova 
feedback. Here, we give a brief outline of the most important 
aspects of the physics mode l ; a m ore complete description 
can be found in lCrain et al. | (|2009l) and references therein. 

• Radiative cooling is computed assuming the optically 
thin limit and ionisation equilibrium, depending on redshift, 
local temperature and density, and chemical composition. 
The IGM is heated uniformly, with a redshift-denenden t 
UV/X-ray background as given bv lHaardt fe Madaul (|200ll ). 
with hydrogen and helium-II being reionised at z ~ 9 and 
z ~ 3.5, respec tively. 

• Following ISchave fc Dalla Vecchial (|2008l ). star forma- 
tion is implemented by imposing an effective equation of 
state for gas particles above a density threshold of nn = 
0.1cm -3 , that makes the Jeans mass independent of the gas 
density. In this regime, gas particles can be transformed into 
star particles at a pressure dependent rate that reproduces a 
local Kennicutt-Schmidt law, with a star formation rate col- 
umn density of E = 1.5 x 10 4 MQyr _1 kpc _2 Eg 4 , where E 9 is 
the gas column density in M0PC -2 . Each star particle is as- 
sumed to represent a sing le stellar populatio n, with an IMF 
based on IChabrierl (|2003l ). As described in IWiersma et all 
120091 ). stellar evolution follows the abundances of 11 metals, 
which are released over time to the surrounding gas particles 
by supernovae type II, type la, and AG B stars . 

• Following iDalla Vecchia fc Schavel (|200Sh . supernova 
feedback is imparted on the surrounding gas particles in 
the form of kinetic energy, with each star in the range of 
6 — 100M© providing 10 jl ergs, sufficient to accelerate an 
average of 4 neighbours to a wind velocity of 600 kms" 1 . 

Whilst some parts of the astrophysical model are still 
subject to considerable uncertainty, both in terms of the 
underlying physics and its numerical implementation, it 
has been shown to reproduce the cosmic star formation 
rate (|Crain et al.l 120091 ; ISchave et al.l |2010|). and the prop- 
erties of individual ga laxies in some detail l|Font et al]|201ll . 
iMcCarthv et alj|2012j ). By comparison, the results presented 
in this paper are quite general. 

The GIMIC simulation does not include AGN. This 
omission is the most likely cause of a too high stellar and 
baryon fraction in the most massive haloes. Because the to- 
tal abundance of these massive haloes is low, the effect on 
the cumulative mass function is negligible for less massive 
haloes. 
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2.2 Analysis 

We compare the populations of haloes and subhaloes in 
the two simulations, as well as individual objects which are 
matched across both. In Section [2.2.11 we describe the iden- 
tification of structures and substructures, and the definition 
of different subsets used in the analysis. In Section 12.2.31 
we describe the procedure and criteria used for matching 
substructures across the two simulations. 



2.2.1 Identification of Substructures 

In each snapshot, structures are identified in a two step pro- 
cess. First, dark matter haloes are found using the Friend- 
of-Friend (FoF) a lgorithm, with a linking length of 0.2 
(|Davis et al.lll985l) . In the GIMIC simulation, star parti- 
cles and gas particles are then attached to the nearest dark 
matter particle, and inherit its FoF association. In a second 
step, self-bound substructures within haloes are ide ntified 
using the SUBFI ND algorithm of Springel et alj (|200lf ) . with 
the extension of iDolag et al.l (|2009l ) to account for baryons. 
Subfind computes the density field within each FoF halo us- 
ing an SPH-interpolation. A set of N > 20 particles (of any 
type) is considered as a potential substructure if it exceeds 
the local smoothed density estimate. It is then subjected to 
unbinding, whereby all particles whose combined kinetic and 
internal energy exceeds their gravitational binding energy to 
the substructure are iteratively removed, until the substruc- 
ture either vanishes, falling below the 20 particle threshold, 
or is identified as a genuine self-bound subhalo. 



2.2.2 Mass Estimators 

Having identified haloes and sub-haloes, their masses can be 
measured in several ways. We define the FoF-mass, Mp p, 
of a halo as the sum of the masses of all particles belonging 
to the FoF-group. In structure formation theory, it is cus- 
tomary to characterise a halo by M200, which describes an 
object with an overdensity of 200, near the threshold for a 
dark matter halo to collapse and virialise. Different defini- 
tions exist for M200 and r2oo, the radius which encloses a 
spherical volume of corresponding overdensity: IVboo.crit is 
defined such that the mean overdensity within r2oo is 200 
x the critical density; M2oo,moan such that the overdensity 
within r2oo is 200 x the mean background density; M top hat 
is the mass in a spherical volume whose overdensity equals 
the value at virialisation in the top-hat collapse model. In 
Fig. [21 we show a comparison of the three mass estimators 
relative to the FoF-mass for all haloes at z = 0. We find that 
Mtophat most closely matches Mp F, and that the spread be- 
tween the three estimators increases with halo mass in the 
GIMIC simulation, from - 20% at 10 9 M Q to ~ 30% at 
10 14 Mq. Importantly for our purpose, below 10 12 Mq the 
ratio between each of the three halo mass estimators and 
the FoF-mass is very similar in both simulations, and much 
smaller than the measured baryonic effect on the halo mass. 
We use Mpop as our halo mass definition throughout the rest 
of the paper, as it is independent of the concentration of the 
halo. The fact that our results are very similar for the dif- 
ferent mass estimators indicates that a measured halo mass 
difference is a true physical effect. 

For subhaloes, the mass Msh is defined as the total mass 



of particles bound to the subhalo. For isolated subhaloes, i.e. 
central subhaloes with no satellites, we find that the subhalo 
mass is typically very close to the halo mass, indicating that 
most particles belonging to the halo are also gravitationally 
bound to it. We have also examined individual subhaloes, 
and compared density profiles of the gravitationally bound 
particles to those of all particles in the vicinity. We find no 
noticeable difference between the two simulations in the ra- 
tio of bound to unbound particles for each subhalo, indicat- 
ing that any subhalo mass difference that we find between 
the two simulations is a genuine physical effect. Examples 
of density profiles, including a comparison of the bound and 
unbound particles, are shown in Fig. [5] 

We discuss the convergence of the simulations in terms 
of the halo- and the subhalo mass functions in Appendix iBl 
and show that the principal effects of the baryon physics 
are independent of the definition of the objects, and are not 
affected by resolution. 



2.2.3 Matching Substructures Across Simulations 

In order to understand how the differences between the 
DMO and the GIMIC simulations arise, we also compare 
the evolution of individual objects. We use the fact that 
dark matter particles of identical position in the initial con- 
dition^] have identical IDs in order to match substructures 
across both simulations. Specifically, at each snapshot and 
for every subhalo in the DMO simulation, we consider its 
5 most bound particles, and require that at least 3 of them 
belong to a single subhalo in the GIMIC simulation. This en- 
sures that each subhalo in the DMO simulation is matched 
to at most one subhalo in the GIMIC simulation. Because 
mergers can occur at different times in the two simulations 
(and in some cases occur only in one and not in the other), 
we also impose a threshold on the mass ratio, matching only 
objects whose mass ratio is less than 4:1. 

As shown in Fig. [3] this method results in > 90% of sub- 
haloes above 1O 9 M0 in the DMO simulation being matched 
to unique counterparts in the GIMIC simulation, including 
> 95% of centrals, and > 80% of satellites. The reduction 
of the matching fraction for present-day satellites compared 
to central subhaloes is largely due to their more chaotic or- 
bits, and the associated increase in the likelihood of diverg- 
ing merger histories of individual objects in the two simula- 
tions. In addition, low-mass satellites are more likely to be 
disrupted in the GIMIC simulation compared to the DMO 
simulation. 

In total, at z = 0, the DMO simulation contains 214,676 
haloes and 254,305 subhaloes; the GIMIC simulation con- 
tains 148,765 haloes and 162,866 subhaloes. Including only 
subhaloes with a total mass above 1O 9 M , the DMO sim- 
ulation has 86,137 subhaloes, 69,699 centrals and 16,438 
satellites, while the GIMIC simulation has 62,265 subhaloes; 
50,077 centrals and 12,188 satellites. Of the central sub- 
haloes in the DMO simulation with masses above 1O 9 M , 
92.6% are matched to subhaloes in GIMIC, 97% of which 
are matched to centrals and 3% to satellites. In the same 



2 A DM particle in the DMO simulation is replaced by a pair of 
DM and gas particles with the same centre of mass in the initial 
conditions of the GIMIC simulation. 
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Figure 4. Projected density distributions of gas (left), stars (centre) and dark matter (right) in the GIMIC simulation, from a cubic 
region of side length I = 11.2 Mpc at z = 0. On large scales, stars and gas follow the dark matter distribution, but on galaxy scales, the 
distributions differ, leading to a difference in the total mass distributions between the GIMIC and the DMO simulations. 



mass range, 80.0% of satellites in the DMO simulation are 
matched to subhaloes in GIMIC, 70% of which are matched 
to satellites and 30% to centrals. 

We also compared the matching success rate for sub- 
haloes in the central parts and near the boundary of the 
high-resolution region (see Section T2.ip . but found no signif- 
icant edge-effects. While the criteria for matching subhaloes 
appear somewhat arbitrary, we note that the matching of 
objects across the two simulations only serves to understand 
the origin of the differences, but that statistical quantities, 
including mass functions and satellite mass functions, are 
independent of the matching rate and criteria. 



3 EFFECT OF BARYONS ON STRUCTURES 

In this Section, we describe the effects that the inclusion of 
baryons have on the formation and evolution of structure 
mainly in statistical terms. We compare the evolution of in- 
dividual objects across the two simulations in Section 13.11 
The effect on the overall abundance of haloes and subhaloes 
is discussed in Section 13.31 Section 13.41 shows the effect on 
the n m ax function, and compares our results to measure- 
ments of the ALFALFA HI survey in the Local Volume. In 
Section 13.51 we focus on the mass functions of individual 
groups and clusters. We highlight the importance of dark 
subhaloes in Section [3.61 which will be relevant to the dis- 
cussion of abundance matching in Section [4] 

From the similarity in the overall mass distribution of 
the dark matter only ("DMO") simulation and its counter- 
part with baryons ("GIMIC") shown in Fig. [1] it is clear 
that the large scale evolution of structures is nearly iden- 
tical. On scales of several Mpc, gravity remains the only 
relevant force even in the hydrodynamical simulation, and 
the gas distribution closely follows that of the dark mat- 
ter. Fig. [J] shows the projected density distribution of gas 
(left), stars (centre) and dark matter (right) in a region of 
11.2 3 Mpc 3 from the GIMIC simulation at z = 0. On scales 
of several hundred kpc and below, it is apparent that the 
gas distribution is much smoother than the dark matter, a 
reflection of the fact that most of the interstellar medium 



of small haloes has been expelled, while the intergalactic 
medium is heated thr ough the UV backgr ound, resulting in 
a large Jeans length jTheuns et alj|2000t ). While the stel- 
lar mass is also highly clustered, it is only significant in the 
most massive objects. 

3.1 Individual Objects 

We use pairs of individual matched subhaloes in the DMO 
and the GIMIC simulation, to study the effect of baryons 
on the formation of haloes and subhaloes. In Fig. we 
show the mass profiles in GIMIC at z = 0, for three typ- 
ical matched central subhaloes of 10 10 , 10 11 and 10 12 Mq. 
For each subhalo, the density in dark matter is shown in 
black, the density in stars in green, and the density in gas 
in blue. As the subhalo mass increases from left to right, 
the baryon fraction also increases. Whereas the total den- 
sity is dominated by dark matter throughout for 10 10 and 
1O 11 M0 subhaloes, 1O 12 M0 subhaloes contain a significant 
mass of stars in the centre. For subhalo masses of 10 10 and 
1O 11 M0, the total mass density in GIMIC is systematically 
lower than the mass density of the corresponding object in 
the DMO simulation, which is plotted as a black dashed line. 
Note that for the 10 12 Mq subhaloes, the central mass profile 
is steeper in GIMIC compared to the DMO simulation as a 
result of cooling and adiabatic contraction; this is not seen 
in the smaller haloes. 

Fig. [6] compares the masses of many individual sub- 
haloes across the two simulations, as a function of their mass 
in the DMO simulation, subdivided into centrals (left panel), 
and satellite (right panel). For the purposes of this compar- 
ison, we only include subhaloes that are of the same type 
in both simulations, i.e. satellites matched to satellites, and 
centrals matched to centrals, using the criteria discussed in 
Section [2.2.31 We show a random subset of subhaloes, nor- 
malised so as to have equal numbers per logarithmic mass 
interval. 

Each subhalo is represented by a set of three points 
that denote the different composition of the subhalo in the 
GIMIC simulation: black circles measure only the dark mat- 
ter component, green circles also include the gas, and red 
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Figure 5. Spherical density profiles of 3 pairs of individual matched, central subhalocs. In the top row, solid lines show the different 
components in the GIMIC simulation: dark matter (black), stars (green), gas (blue), and total mass density (red), each including only 
particles bound to the subhalo. The dashed line shows the density of the same object in the DMO simulation, and vertical lines indicate 
the radius of the outermost bound particle. The material within ~ 1O 12 M0 subhaloes has been rearranged appreciably in the GIMIC 
simulation relative to the DMO simulation, with baryons accumulating in the centre. The difference between the solid and dashed black 
lines in the inner part shows that the dark matter itself is also more concentrated in GIMIC, but the total subhalo mass does not change. 
By contrast, ~ 10 10 Mq and ~ 10 11 Mq subhaloes have lost mass in the GIMIC case. Here, the total density is dominated by dark matter 
throughout, and lower than in the DMO simulation. In the bottom row, solid lines denote the density profiles computed using bound 
particles, dotted lines show unbound particles in the vicinity (not counted towards the subhalo mass), and dashed lines show the sum of 
the densities of bound and unbound particles. In all cases, the density of unbound particles is low, and there is no significant difference 
between the two simulations in the ratio of bound to unbound particles. 




Figure 6. Ratio of subhalo mass in GIMIC relative to the DMO simulation for individual, matched subhaloes at z = 0, as a function 
of DMO mass, for central subhaloes (left) and satellite subhaloes (right). Every subhalo is represented by a set of three circles, with the 
y- values indicating its mass components in the GIMIC simulation: only the dark matter (black), dark matter and gas (green), and total 
mass including dark matter, gas and stars (red). The solid lines of corresponding colour indicate the median value within each bin, from 
which individual subhaloes arc drawn at random, so as to show the same number of objects per logarithmic mass interval. The black 
dashed and dotted lines denote a 1:1 ratio, and a ratio of CI dm I {^DM + ^b)> respectively. For lower mass subhaloes, the circles lie closer 
together, indicating a decrease in baryon fraction with decreasing subhalo mass, reflected also in the convergence of the three curves. 
The decline of the red curve indicates a decrease in total subhalo mass in GIMIC relative to the DMO simulation. For satellites, there 
is a similar trend of decreasing baryon fraction and mass ratio for lower mass objects, although the scatter among individual subhaloes 
is much greater. 
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circles represent the total mass, consisting of dark matter, 
gas and stars. 

Overplotted on both panels are the median mass ratio of 
the corresponding components; black lines for dark matter, 
green also including gas, and red including dark matter, gas 
and stars. As the baryon fraction decreases with decreasing 
subhalo mass, circles which are well separated at high mass 
become closer, and for lower mass subhaloes, they mostly 
overlap. This trend is reflected in the convergence of the 
median mass ratios of the three components. At the high- 
mass end, on average, subhaloes reach the same total mass 
in the GIMIC simulation as in the DMO simulation because 
their baryon fraction is close to the universal value. The total 
mass (red line) approaches a 1:1 ratio, while the dark matter 
component (black line) approaches the value expected from 
subtracting the universal baryon fraction. At lower masses, 
the median total mass in the GIMIC simulation is dominated 
by the dark matter component. It is worth noting that the 
mass of low-mass subhaloes in the GIMIC simulation is not 
only below the mass in the DMO simulation, but even be- 
low the value expected from a simple loss of baryons. This 
points to the fact that objects which lose mass during an 
early stage in their formation acquire a shallower potential, 
which also diminishes their subsequent accretion. Due to the 
requirement for subhaloes to be self-bound, a subhalo whose 
gravitational binding energy is reduced due to outflows may 
also lose additional particles that are no longer bound. How- 
ever, as we show in Fig. [9] we find a very similar reduction 
in the masses of FoF-haloes, indicating that any difference 
in the attribution of masses by SUBFIND is only a secondary 
effect. 

Comparing the central subhaloes in the left panel of 
Fig. [5] to the satellite subhaloes in the right panel, it can 
be seen that the scatter increases substantially for satellites, 
and also increases with decreasing mass. As we will show in 
Fig. 1121 the overall gas fraction is lower for satellites than 
for centrals as a result of stripping. 

The main cause for the increased scatter among satel- 
lites lies in their more chaotic orbits, meaning that matched 
pairs are more likely to follow divergent evolutionary paths 
in the two simulations. This effect is illustrated in Fig. [JJ 
which shows the projected displacement between pairs of 
satellite subhaloes of ~ 1O 1O M0 in the GIMIC simulation 
relative to the DMO simulation. Tracing the evolution since 
2 = 1.5, it can be seen that satellite pairs whose final mass 
ratios deviate from the median by more than 2a (top panel) 
are much more likely to evolve along divergent paths com- 
pared to those whose final mass ratios are within 10% of the 
median (bottom panel). 

Qualitatively, the overall reduction in subhalo mass due 
to baryonic effects at z — is similar for centrals and satel- 
lites. However, despite the considerable scatter, at fixed sub- 
halo mass, the average mass loss in the GIMIC simulation 
appears to be somewhat larger for central subhaloes than for 
satellites, as is also illustrated in Fig.[8]and Fig. lAll . The off- 
set between the relations seen for satellites and for centrals 
may be partly attributed to the fact that while centrals only 
experience mass loss due to outflows in the GIMIC simula- 
tion, the masses of satellites are further reduced by stripping 
in both the DMO and the GIMIC simulations. 

It is worth noting that, despite the fact that every halo 
contains at most one central subhalo but potentially many 




Figure 7. Relative displacement between matched pairs of satel- 
lite subhaloes with final masses of ~ 10 10 Mq, for pairs with mass 
ratios that deviate from the median by more than 2<r (top panel) , 
or by less than 10% (bottom panel). Circles show the present dis- 
placement, while traces show the evolution from z = 1.5, coloured 
according to the mass ratio: black for mass ratios close to the me- 
dian, red (blue) if the mass in GIMIC is higher (lower) than the 
median at the time. While most pairs of subhaloes in the bottom 
panel evolve along similar paths, those in the top panel show more 
divergent histories, linking the increased scatter among satellites 
in Fig. [6] to their more chaotic orbits. 



satellites, over 80% of all subhaloes in the simulation volume 
are central, as was shown in Fig. [3] This is easily understood 
when recognising that within each halo, by definition, the 
central subhalo is more massive than the satellites, so at 
fixed subhalo mass, satellites belong to larger haloes than 
centrals. As a consequence of the steep halo mass function 
(see Fig. [9]), the total subhalo population closely resembles 
that of centrals, with some additional scatter introduced by 
satellites. 
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Figure 8. Ratio of subhalo masses between the GIMIC and DMO 
simulations for matched subhaloes at z = compared to analytic 
fits. The blue and red points show the results for centrals and 
satellites, respectively, while the error bars show the estimate of 
the median ratio and its statistical error, for centrals (blue), satel- 
lites (red) and all subhaloes (black). The solid and dashed lines 
show the best fit to the median for all subhaloes, assuming 4 free 
parameters (dashed line, Eq. [TJ, or 3 free parameters, with the 
upper limit fixed at b = 1. (solid line, Eq.[2} 



3.2 Analytic Correction 

Noting from Fig.[6]that the change in mass ratio between the 
GIMIC and the DMO simulations saturate both at the high- 
mass and at the low-mass limit, we parametrise the mean 
ratio of the subhalo mass in the GIMIC simulation, Mgimic, 
to that in the DMO simulation, Mdmo , in the following way: 



Mgimic _ b a/b + (Mdmo /M t ) 



Mdmo 1 + (Mdmo /Mt ) w 

With w > 0, Eq. [T]has two horizontal asymptotes, 

/ Mgimic 1 



(1) 



lim 

M DMo<< M t V Mdmo 



and 



lim 

M D MO»Mt 



/ Mgimic \ 
t V Mdmo / 



= b. 



The mass ratio takes an intermediate value of (a + 6)/2 
at Mdmo = Mt, i.e. Mt sets the scale for the transition 
between the two limits, while the width of the transition 
is parametrised by the value of w. Assigning equal weight 
to the median ratio in each logarithmic mass interval, the 
best fit to all subhaloes at z = yields values of a — 0.69, 
b = 0.98, M t = 1O 11<3 M and w = 0.79. It is shown by the 
dashed line in Fig. [8] 

The asymptotic behaviour in the low-mass limit signi- 
fies an almost complete loss of baryons at very early times. 
The fact that the value of b at z = is very close to unity 
indicates that the physical processes in GIMIC lead to no 
significant change in mass at the high mass end. By setting 
6=1, Eq. [1] simplifies to: 



Mgimic _ a + (MoMo/Mt) 1 



Mdmo 1 + (M DM o/M t ) 1 



(2) 



a = 0.65, M t = lO n - 4 M and w = 0.51. The solid line in 
Fig. [8] shows the fit to Eq. [2] and can be compared to the 
dashed line obtained from Eq. [T] 

Because the additional free parameter of Eq. [T] does 
not change the results significantly, we adopt the simpler 
Eq.[2]for the remainder of this work. It is possible that other 
physical processes, particularly the effect of AGN, could lead 
to a significant effect also at the high mass end, in which case 
the more general form of Eq. [T] would be more appropriate. 
Separate fits for satellites and centrals, as well as for different 
redshifts, are discussed in Appendix lAl 



3.3 Global Mass Functions 

Mass- or multiplicity functions of haloes or subhaloes con- 
stitute one of the most fundamental measures of structure 
formation. Because dark matter haloes (or more precisely: 
subhaloes), are believed to be the locations of galaxy for- 
mation, predicting their abundance also allows us to relate 
observations of galaxies to theories about the underlying 
mass distribution. We return to this point in Section [4] 

In Fig. [5] we show the cumulative mass functions of 
haloes (left panel) and subhaloes (right panel), for the DMO 
simulation (shown in black) and the GIMIC simulation 
(red). In both cases, the mass includes the total mass; i.e. 
the mass in dark matter in the DMO simulation, and the 
combined mass of dark matter, gas, and stars in the GIMIC 
simulation. The cumulative mass functions are nearly iden- 
tical for the most massive objects, but below ~ 10 12 Mq, 
the abundance in the GIMIC simulation is significantly re- 
duced relative to the DMO simulation, with the difference 
increasing at lower masses. 

It should be noted that this statistical result is not re- 
liant on the matching of haloes or subhaloes and, as we show 
in Appendix [B] it is not affected by the resolution of our 
simulations. Above the resolution threshold, the decrease in 
numbers is only partly due to the erasure of structures, but 
mostly due to each object having a lower mass in the GIMIC 
simulation. The shape of the mass functions then results in 
this "horizontal" shift in mass also manifesting itself as a 
"vertical" shift in abundance. 



3.4 Velocity Functions 

The circular velocity, v c (r), at a given radius, r, is related 
to the enclosed mass m(< r) via v c (r) = U Gm(< r) jr. 
Thus, measurements of the velocity profiles of HI gas in 
galaxies can provide a measurement of the underlying mass 
distribution. In Fig. [10] we compare the differential veloc- 
ity function of subhaloes in the DMO simulation to results 
from the GIMIC simulation, and to measurement s from the 
Arec i bo Le gacy Fast ALFA (ALFALFA) survey (|Giovanellil 
120051 . l2007h. Following constra ined simulations of the Local 



with horizontal asymptotes at a and 1, and an intermediate 
value of (a + l)/2. A fit to this equation yields values of 



Volume by I Zavala et al.l(|2009l ), we have normalised the sim- 
ulated velocity functions to account for the overdensity in 
the survey volume, which is dominated by the Local Super- 
cluster, and whose velocity function exceeds the universal 
value. A detailed des c riptio n of the normalisation can be 
found in lZavala et al.l l|2009l ). but we note that it has little 
effect on scales where the DMO results significantly exceed 
the observations. 
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Figure 9. Comparison of the cumulative halo mass functions (left column), and subhalo mass functions (right column) of the GIMIC 
and the DMO simulations. The top row shows the cumulative mass functions at redshift 2 = 0, with the DMO results shown in black, 
and the GIMIC results overplotted in red. The bottom row shows the ratio of the cumulative abundances at z = 0, with the shaded area 
denoting the statistical uncertainty. The difference in abundance depends strongly on mass: Above ~ 1O 12 M0, both the subhalo and 
halo abundances have ratios close to 1, but they drop to ~ 0.75 from ~ 10 n MQ and below. 



When compared to the ALFALFA observations, it can 
be seen that the velocity function of the DMO simula- 
tion (black) over-predicts the abundance of subhaloes of 
35kms _1 by a factor of ~ 10. This d iscrep - 



ancy, which has been point ed ou t by | Zavala et al 



iTruiillo-Gomez et al.l (|201ll ) and IPapastergis et al 



2009), 



20111), 



persists when the observations are compared to the veloc- 
ity function computed using all subhaloes in GIMIC (red). 
It has been considered a challenge to the cold dark mat- 
ter (CDM) paradigm, and evidence in favour of warm dark 
matter (WDM). However, allowing for the detection thresh- 
old oftJ^jsurvey for th e sample of nearby galaxies analysed 
bv lZavala et all J20091 I, Mm > lO 7 h~ 1 M , and recomput- 
ing the GIMIC velocity function only for subhaloes with 
M g > lO 7 h _1 M0, we find that the resulting velocity func- 
tion (green symbols) agrees to within a factor of ~ 3 with the 
ALFALFA observations, while also giving a better fit to the 
shape of the t) max function. Considering that the total gas 
mass, M g , is only an upper limit for Mm, and that we have 
not included other constraints such as the inclination limit 
of 30 for galaxies in the ALFALFA survey, the remaining 
difference may be yet be explained within the CDM frame- 
work. 



3.5 Satellite Mass Functions 

Local Group dwarf galaxies and Milky Way satellites consti- 
tute important test cases for CDM theory: they are among 



the smallest and most dark-matter dominated galaxies, and 
their proximity enables detailed studies of their abundance 
and their internal structure. While the apparent discrep- 
ancy between the predicted satellite mass function and 
the observed satell i te lum inosity function (e.g. lKlypin et al] 
1 19991 : iMoore et aljfl999h may be explained through ineffi- 
cient star formation due to astrophysical mechanisms includ- 
ing photo-ionisatio n and feedback (e .g. iBullock et al. I l200d : 
iBenson et alj [20021 : ISomervilIell2002r i. the reported discrep- 
ancies between the predicted and measured satellite mass 
functions are more challenging: ACDM s imulations such as 
the Aquarius runs (|Springel et alj |200S( ) produce a higher 
number of satellites with circul ar velocities above 50 kms" 1 
than is allowed by obser vations (|Bovlan-Kolchin et alj|2()lll . 
120121 : IParrv et alj|2012r ). It has also been pointed out that 
the Milky Way appears special in ha ving two sa t ellites 
as bright as the Magellanic Cloud s (|Busha et alj 1201 ll : 
iTollerud et~ai1l201ll ; IGuo et alj|201ll ). 



In Fig. 111! we compare the satellite mass functions of 
groups of different mass in the DMO simulation (black lines) 
and the GIMIC simulation (red lines). In the top panel, in 
both simulations, we have selected groups with total mass 
in the range of 1 — 2.5 x 1O 12 M0, comp atible with current 
mass estimates of the Milky Way (e.g. iLi fc White] 120081 : 
IXue et alj [20081 : iGnedin et al]|2010h . Each thin line shows 
the satellite mass function of an individual group, while 
the two thick lines represent the average cumulative sub- 
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Figure 10. Differential velocity function of subhaloes in the 
DMO simulation and in the GIMIC simulation, compared to the 
HI velocity function measured in the ALFALFA survey for a de- 
tection limit of Mm = 10 7 h _1 M Q . The results of the DMO simu- 
lation are shown in black; for GIMIC, red shows the velocity func- 
tion computed using all subhaloes, and green only includes those 
subhaloes with M g > lO 7 h _1 M0. The simulated mass functions 
are corrected for the overdensity of the survey volume, their width 
denote the Poisson errors within each bin. The vertical lines at 
20 kms - 1 and 35 kms - 1 indicate the resolution limit of our sim- 
ulations and the completeness limit of the ALFALFA survey, re- 
spectively. While the DMO simulation significantly over-predicts 
the abundance of subhaloes with v ma x = 35kms -1 , accounting 
for baryonic effects and adopting an appropriate detection thresh- 
old for the ALFALFA survey significantly reduces the reported 
discrepancies between CDM and the measured velocity function, 
down to the survey limit of 35 kms" 1 . 



halo abundance over all groups in the mass range in both 
simulations. 

We find that, on average, the MW-mass groups contain 
fewer satellites in GIMIC than in the DMO simulation. The 
reduction in the number of satellites in haloes of a fixed mass 
can be understood because the mass-dependent baryonic ef- 
fects are stronger in the lower mass satellites than in their 
more massive hosts. 

In line with previous results bv ldi Cintio et all (|201ll >. 
the average baryon effects we see in the GIMIC simula- 
tion alone are not enough to sufficiently reduce the num- 
ber of massive satelli t es in order to explain the results of 
iBovlan-Kolchin et alj l|201ll ). However, in both the GIMIC 
and the DMO simulation, we also find a scatter of more 
than one order of magnitude in the satellite abundances 
of individual groups. In particular, among the sample of 
~ 50 Milky- Way mass groups, we find several examples with 
satellite numbers almost an order of magnitud e below the 
median. Confirming the results of IWana l|2012T l. we there- 
fore find a considerable chance for a low number of massive 
satellites in a Milky- Way mass halo, to which baryon effects 
contribute further. 
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Figure 11. Cumulative satellite mass functions for ~ 50 individ- 
ual, Milky- Way sized (10 12 M Q < M 200 < 2.5 X 10 12 M Q ) haloes 
(top), and for the 20 most massive clusters (bottom), in the DMO 
simulation (black) and the GIMIC simulation (red), at z = 0. 
Thin lines show the mass functions of individual groups, thick 
lines show the average abundance for the groups in either sim- 
ulation. For clusters and for Milky Way groups, the number of 
members is reduced in the GIMIC simulation, as expected from 
the overall reduction in subhalo abundance. It is also worth not- 
ing that there is more than an order of magnitude variation in 
the number of satellites amongst individual groups. 



In the bottom panel of Fig. 1111 we repeat the same 
analysis for the satellite mass functions of the 20 most mas- 
sive groups in both simulations; clusters with masses in the 
range ~ 10 13 — 10 14 Mq. As expected from the global mass 
functions shown in Fig. [9] at the high-mass end, the aver- 
age satellite mass functions of clusters show no systematic 
differenc e between the two simulations, a result noted previ- 
ously bv lDolag et~aj] i|2009f ). However, at lower masses, clus- 
ters also have fewer satellites in GIMIC than in the DMO 
simulation. Compared to the satellite mass function of MW- 
sized groups shown in the top panel of Fig. 1111 the satellite 
mass functions of individual clusters show less scatter, even 
though the mass range of clusters shown is larger than that 
of MW-sized groups. 
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Figure 12. Fraction of subhalocs with stars (blue), with gas (red), or both (black), at z = 0. Solid lines show the results for centrals (left 
panel) and satellites (right panel); dashed lines show the fraction of all subhaloes for reference and are identical in the two panels. For 
both centrals and satellites, the fraction of subhaloes with stars decreases for subhalo masses below ~2x 10 10 Mq, while the fraction 
of subhaloes with gas starts to decrease below ~ 3 X 10 1:l Mq. The main environmental dependence is in the gas content; while the 
fraction of central galaxies with gas remains above ~ 50% down to 10 9 Mq, it drops below ~ 25% among satellites with a total mass of 
~ 2 X 10 10 Mq. The decrease in the fraction of subhaloes without stars is less steep for satellites than for centrals; this is because stars 
form in more bound objects, which are more likely to survive as satellites. 



3.6 Dark Subhaloes 

At z = 0, the GIMIC simulation contains a substantial num- 
ber of low-mass subhaloes which are completely free of stars 
or gas, i.e. that they do not contain a single (bound) star 
or gas particle of mass 1.98 x 1O 6 M0. Dark subhaloes are 
predicted to exist in ACDM, as they may explain the ratio 
of the predicted number of substructures to the observed 
number of dwarf galaxies I Moore et al.l (|l999l ) , and because 
astrophysical processes including reionisation prevent cool- 
ing and star formation in the smallest systems. Likewise, 
very low-mass galaxies like dwarf spheroidals contain stars 
but no detectable gas, which is understood to have been 
lost either through supernova feedback, reionisation (if stars 
formed earlier), stripping, or a combination of these effects. 

In Fig. [12] we show the fraction of low-mass subhaloes 
in the GIMIC simulation containing some amount of stars 
(blue lines), gas (red lines), or both (black lines), as a func- 
tion of total subhalo mass. For both centrals (left panel) 
and satellites (right panel), we find that the fraction of dark 
subhaloes increases with decreasing mass below ~ 10 10 M Q 
(i.e. ~ 10 3 particles). Intriguingly, the fraction of subhaloes 
with stars decreases less steeply with mass among satel- 
lites compared to centrals; while ~ 50% of satellites with 
a total mass of 1O 9 M0 contain stars, the fraction is below 
10% for centrals of the same mass at z = 0. This difference 
can be understood from the effective threshold in binding 
energy required for star formation, apparent from Fig. 1131 
subhaloes that are more concentrated, and therefore have a 
higher maximum circular velocity v max , are both more likely 
to cool gas efficiently to allow star formation, and also less 
likely to be destroyed by stripping. As can be seen in Fig. 1131 
a Vmax threshold of ~ 30 kms , rather than a threshold 
in total mass, separates subhaloes that contain stars from 
those that do not. In preferentially destroying less bound 
objects, tidal disruption leads to a higher fraction of sub- 



haloes with stars among the surviving satellites relative to 
the fraction among centrals, which suffer weaker tidal forces. 
Our r esult is in agreement with simulations of lCrain et al.l 
|2007l ). who found that photo- heating can prevent the col- 
lapse of baryons in systems with circular velocities below 
35 kms ~ (at z = 0), and wit h the high-resolution simula- 
tions of lOkamoto et al I l|2008l ). who found that haloes with 
a circular velocity of 25 kms -1 can lose half of their baryons 
due to photo-heating, and are unable to accrete and cool gas 
efficiently. 

The fraction of subhaloes with gas decreases from ~ 
1O 11 M0 downwards, and decreases much faster for satellites 
than for centrals. While stars in satellites are concentrated 
towards the centre, gas is easily re moved from satellites b y 
tidal and ram pressure stripping (e.g. lMcCarthv et alj20Cllil ). 
At 10 10 Mq, where most subhaloes contain stars, the fraction 
with gas is ~ 55% among centrals, and only ~ 20% among 
satellites. 

We also find that the fraction of subhaloes without gas, 
and the fraction without stars, are correlated, but not pro- 
portional. The black lines in Fig. [15] show the fraction of 
subhaloes that contain both stars and gas. It can be seen 
that central galaxies with total mass below ~3x 1O 9 M 
typically only contain either gas or stars: subhaloes that 
undergo star formation lose their gas due to feedback, while 
those that never form stars due to a lack of cooling, can 
typically retain their gas. 

As subhaloes without stars are not directly observable, 
their presence has important consequences for abundance 
matching, that will be addressed in Section [4] Of course, 
the finite particle mass of m g = 1.98 x 10 6 Mq sets a natural 
lower limit to the baryonic mass detectable in our simula- 
tion. However, the occurrence of dark subhaloes at relatively 
high subhalo masses, and the dependence on v max and envi- 
ronment, rather than on subhalo mass, indicate the physical 
origin of their existence. As we will discuss in Appendix ITU 
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Figure 13. Distribution of maximum circular velocities, v max , of 
a random sample of low-mass subhalocs in the GIMIC simulation, 
as a function of subhalo mass. Subhaloes that contain stars are 
shown in red, those without stars are shown in black At masses 
below ~ 10 10 , it is apparent that a v max of ~ 30 kms" 1 splits the 
two populations: subhaloes with v ma x > 30 kms -1 have typically 
been able to form stars, while most subhaloes with lower v ma x 
have not. 



while the gas-free subhaloes are well converged, predictions 
for the exact fraction of dark subhaloes are still affected by 
resolution in the GIMIC simulation. 



4 ABUNDANCE MATCHING 



Abundance match i ng (lYang et al 



l200rj ; iMoster et al.|[2010l ; IGuo et all 



2001; IVale fc Ostrikeil 
20101 ). is a simple but 



powerful method to statistically link galaxies to dark mat- 
ter haloes, without making detailed assumptions about the 
physics of galaxy formation. By assuming a monotonic re- 
lationship (with a limited amount of scatter) between an 
observed property such as stellar mass and a simulated prop- 
erty such as total mass, and equating their cumulative abun- 
dances, it allows a determination of quantities like the stel- 
lar mass-total mass ratio as a function of mass. In integral 
form, the relation (for stellar mass) can be expressed, in the 
absence of scatter, as follows: 



Nh(m)drn = 



N+(m)dm, 



(3) 



where Nh(m) and N+{m) are the subhalo and stellar mass 
functions. The upper limits of the subhalo and stellar mass 
ranges, Mh, ma x and M* >max , are fixed, while the respective 
lower limits, Mh, m in and M*, m i n , are chosen such that the 
relation holds over the entire interval. Equivalently, abun- 
dance matching corresponds to assigning the most massive 
observed galaxy to the most massive (sub)halo, and step- 
ping down towards the minimum galaxy mass or subhalo 
mass that can still be matched, which determine Mh, m i n and 
M*. min . 




Figure 14. Mass of subhaloes in the GIMIC simulation and in 
the DMO simulation, at equal cumulative subhalo abundance. The 
blue line includes all subhaloes, while the red line assumes the re- 
duced mass function, with dark subhalocs removed. The dashed 
line provides a reference for a 1:1 abundance ratio. Below 10 12 Mq , 
the mass attributed to subhaloes of equal abundance is consid- 
erably lower in the GIMIC simulation compared to the DMO 
simulation, and the occurrence of dark subhaloes becomes the 
dominant factor below ~4x 10 9 Mq. 



4.1 Abundance matching with baryons 

The effect of baryons on the abundance of subhaloes shown 
in Fig. [9] has consequences for the results of abundance 
matching. Furthermore, as described in Section 13.61 an in- 
creasing fraction of low-mass subhaloes does not host any 
stars. Because these would not be observable, we propose 
to modify Eq. [3] further, by multiplying the (sub)halo mass 
function by a completeness term ^ f*(jn) ^ 1, to give the 
reduced subhalo mass function (RMF) of observable objects: 



f i ,(m)Nh(m)dm — 



N i ,(m)dm, 



(4) 



Unlike a simple scatter term, dropping the assumption that 
every simulated dark matter subhalo with m > Mh.min hosts 
a galaxy has the effect of extending the domain to lower mass 
subhaloes, and assigning higher stellar masses to subhaloes 
above the new limit. 

For the purpose of abundance matching, the reduced 
subhalo mass function differs depending on the observable: 
for a stellar mass survey, the RMF includes o nly subhaloes 
that contain stars, while for an HI survey (e.g. IZavala et all 
l2009h . it includes those galaxies that contain gas, or both 
stars and gas. 

From the ratio of the subhalo mass functions, we com- 
pute the subhalo mass in GIMIC with the same cumulative 
abundance as a subhalo of a given mass in the DMO sim- 
ulation. Fig. [14] shows the relations of equating the DMO 
subhalo mass function to the total subhalo mass function in 
GIMIC (blue line), and to the reduced subhalo mass func- 
tion in GIMIC (red line), where subhaloes without stars have 
been removed. These relations can be applied directly to 
abundance matching results based on a dark matter only 
simulation, to derive stellar-to-total mass ratios in a uni- 
verse that contains baryons. 
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Figure 15. Stellar mass as a function of subhalo mass (left panel), and ratio of stellar mass to subhalo mass (right panel). Black lines 
are adapted from Guo ct al. (2010), using the halo mass function of the Millennium I and II dark matter only simulations. Blue lines 
are the result of taking into account the baryonic effects on the subhalo mass function measured in GIMIC, as shown in Fig. QJ] Red 
lines are obtained using the reduced subhalo mass function (RMF) of GIMIC, with dark subhaloes removed. Solid lines correspond to 
a faint-end slope of the stellar mass function o f a = —1.15 (Li & Whi te 2009), while dashed lines assume a = —1.58. For comparison, 
we show stellar masses (Misgcld & Hilkcr 2011) and dynamical masses (Pcfiarrubia ct al. 2008) of Milky Way satellites. The inclusion of 
baryons doubles the inferred stellar-to-total mass ratio for subhaloes below 1O 11 M0, while the effect of the RMF further increases the 
inferred ratios by up to an order of magnitude for subhaloes in the range 10 9 — 10 10 Mq. Combined with a steep faint-end slope for the 
stellar mass function, the reported discrepancy between abundance matching and measurements of individual galaxies is alleviated. The 
RMF also results in a finite minimum stellar-to-total mass ratio in subhaloes below ~ IO^Mq. 



Appl ying our re s ults t o the abundance matching per- 
formed bv lGuo et all ()2010l ). we note that the GIMIC simu- 
lation uses the same cosmolo gical param eters as the Millen- 
nium simulations studied by Guo et all . However, while the 
total abundance of (sub)haloes depends on the slope and the 
norm alisation of the power spectrum (e.g. IVale fc Ostrikerl 
120061 ). within the CDM paradigm, the effect of baryons 
should be similar. 

In line with other authors, iGuo et al] have computed 
the subhalo mass function as a combination of the present- 
day mass for centrals, and the mass before infall for satel- 
lites. In adopting their result for the DMO case which forms 
the baseline of all of our abundance matching results, we 
use the same prescription, but we consider the ratio of the 
abundance at 2 — when we apply the baryon effects. As 
shown in Fig. IAU the average mass ratio shows little evolu- 
tion between z — 1 and 2 = 0, suggesting that the difference 
between applying the correction at infall or at 2 = would 
be negligible. 

In Fig. 1151 we show the effect of baryons on the results 
of abundance matching, in terms of the inferred stellar mass 
(left panel), and the stellar-to-total mass ratio (right panel). 
For reference, black lines show the DMO case; th e solid line 
denotes the original results of IGuo et al.l (|2010l ). while the 
dashed line reflects a change in the fai nt-end slope of the 
stellar mass fun ction from a = — 1.15 (|Li fc White! 120091 ) 
to a = -1.58 (|Baldrv et alj|2008h . The set of blue lines 
shows the corresponding relations, but taking into account 
the baryonic effects on the subhalo mass function obtained 
from the GIMIC simulation. The set of red lines uses the 
reduced subhalo mass function (RMF), which includes bary- 
onic effects, and additionally has the dark subhaloes removed 
from the subhalo mass function. 



The effect of baryons on the attributed stellar mass 
strongly depends on the subhalo mass. It is insignificant 
at subhalo masses above 10 12 Mq, but appreciable at the 
faint end. Purely due to the inclusion of baryons, the stellar 
mass-total mass ratio increases by a factor of ~ 1.5 for sub- 
haloes of 1O 11 M0, and by a factor of ~ 2.2 for subhaloes of 
10 10 Mq and below. If the RMF is considered, stellar masses 
assigned to subhaloes below 1O 1O M0 increase further, and 
the mean stellar-to-total mass ratio reaches a finite mini- 
mum at 1 : 10 3 - 1 
stellar mass function. 



10 , depending on the slope of the 



The apparent discrepancy in the stellar-to-total mass 
ratios inferred from abundance matc hing, and state-of- 
the-art simulations of dwarf galaxies (| Sawala et al.l 1201 ll . 
and references therei n), as well as kinem atics of individ- 
ual dwarf galaxies l|Ferrero et al.1 I2012T ). have been de- 
scribed as a challenge to the CDM paradigm. Included in 
the left panel of Fig. fT5] are observational estimates for 
7 cla ssical dwarf spheroidals , with stellar masses adopted 
fromlMisgeld fc Hilkerl <l201ll) . and virial masses derived by 
iPenarrubia et al.1 (|2008l ). 

It is worth pointing out that at such low masses, a di- 
rect comparison to observations should be taken with a grain 
of salt: mass measurements rely on assumptions about the 
shape and anisotropy of the dark matter halo; the halo mass 
is measured at present rather than at infall, and the limited 
sample may be b iased towards the objects with higher stel- 
lar mass (indeed. I Wolf et ail (|2010l ) show that despite stellar 
masses extending to ~ 1O 4 M0, most Milky Way satellites 
are compatible with a formation in ~ 1O 9 M0 haloes). Nev- 
ertheless, it can be seen that when the CDM subhalo mass 
function is corrected to take into account baryonic effects, 



© 2012 RAS, MNRAS 000.1111191 



15 



the stellar masses inferred from abundance matching become 
con sistent with observat i ons of individual dwarf galaxies. 

iTikhonov fc Klypinl l|2009h have likewise argued that 
the velocity function of voids in the Local Volume presents 
a challenge to CDM. The observed number of voids is so 
high that only haloes down to masses of 6 — 8 x 1O 9 M0 or 
1 - 2 x 10 10 M Q (for values of cr 8 of 0.75 and 0.9, respectively) 
can be populated by galaxies brighter than Mb = —12, but 
individual dwarf galaxies are observed with HI rotational 
velocities below iw x ~ 25 kms -1 indicative of halo masses 
closer to 1O 9 M0. Our results suggest that this discrep- 
ancy may also be resolved by baryonic processes. As shown 
in Fig . 1141 the mass limits derived by ITikhonov fc Klypinl 
(j200gh based on the abundance of haloes in dark matter only 
simulations would have to be reduced in a universe where 
baryons affect the masses of low mass haloes. In addition, 
the appearance of dark subhaloes discussed in Section 13.61 
instead of a sharp threshold in (sub)halo mass can reconcile 
the observation of a population of galaxies in lower mass 
haloes with theoretical predictions. 



5 DISCUSSION 

We have compared a set of high-resolution cosmological N- 
body simulations with and without baryons to study the 
effect that baryonic processes have on structure formation 
in a ACDM universe. While both simulations agree well on 
large scales, objects below ~ 10 12 Mq have systematically 
lower masses in the GIMIC simulation that includes baryons 
compared to its "Dark Matter Only" (DMO) counterpart. 
Consequently, the cumulative abundance of haloes and sub- 
haloes at every mass below ~ 10 1 Mq is reduced in the 
GIMIC simulation. Given that the Universe includes bary- 
onic processes such as gas pressure, cooling, reionisation, 
star formation and supernova feedback, and assuming that 
they are represented at least approximately in the GIMIC 
simulation, it follows that DMO simulations over-predict the 
true abundance of substructures, by ~ 10% at lO n ' 5 M , 
~ 20% at lO n M , and ~ 30% at 1O 1O M . 

The difference in halo and subhalo abundance has con- 
sequences for all analyses that take their mass functions as 
a starting point. For example, the stellar-to-total mass ratio 
of dwarf galaxies with stellar masses of ~ 10 5 - 1O 9 M in- 
ferred from subhalo abundance matching doubles after tak- 
ing into account the lower mass of subhaloes, resulting from 
baryonic effects. If, in addition, the reduced mass function 
accounts for the fact that not all subhaloes contain stars, 
stellar-to-total mass ratios increase further, particularly for 
galaxies with stellar masses below ~ 1O 6 M , i.e. the classical 
dwarf spheroidals. It appears that both effects are needed in 
order to reconcile abundance matching results with kine- 
matic estimates of mass-to-light ratios of dwarf galaxies 
(iFerrero et al.ll2012ri. or current high-resolution simulations 



(e.g. ISawala et al.ll20li"h 



We conclude that the discrepancy between dwarf galaxy 
abundances and kinematics, and the results of cold dark 
matter (CDM) simulations paired with abundance match- 
ing, may be due to the inability of dark matter only sim- 
ulations to capture the physical effects relevant on small 
scales, aggravated by the assumption that a visible galaxy 
forms inside every single dark matter (sub)halo. Likewise, 



taking into account baryonic effects and excluding sub- 
haloes that fall below the observational detection threshold 
significantly reduces the reported discrepancy between the 
ACDM paradigm, and the velocity function measured in the 
ALFA LFA HI survey (|Zavala et alj|2009l ; IPapastergis et all 

[2011]). 

For haloes more massive than ~ 1O 12 M , our results 
corroborate tho se of previous works (e.g. iDolag et al. I l2009l ; 
ICui et al.ll20l3 ), who concluded that the net effect of gas 
loss via stripping and feedback is largely counterbalanced 
by the increased concentration of the stellar component. We 
cannot exclude the poss ibility, however, that as suggested by 
Ivan Daalen et ail |201ll ). efficient AGN fee dback would also 
chang e the results on these scales; indeed iMcCarthv et al.l 
|2012l ) have reported overcooling in galaxies with stellar 
masses above 1O 11,5 M in the GIMIC simulation. While 
there are also differences in the physics model, most no- 
tably the inclusion of photo-ionisation in the GIMIC simu- 
lation that suppresses star formation in small objects, the 
most significant difference compared to previous works is 
the increased resolution, which has allowed us to study ob- 
jects more than two orders of magnitude lower in mass. As 
we show, the baryonic results are strongly mass dependent. 
At low masses, re-ionisation, feedback and stripping out- 
weigh adiabatic contraction, leading to a net mass loss. The 
GIMIC results also confirm the results of high-resolution 
resimulations of ind ividual dark matter haloes with baryons 
|Sawala et al.l l201lh . which showed a similar reduction in 
mass due to efficient outflows from ~ 1O 1O M haloes. 

Quantitatively, our results will likely be refined by fu- 
ture simulations at still higher resolution and with more 
complete physics models. In particular, our result for the 
fraction of dark subhaloes straddles the resolution limit of 
our simulations. However, the requirement of a reduced sub- 
halo mass function to match the faint end of observable 
galaxies, of taking into account the gas fraction when com- 
paring to HI surveys, and of including baryon physics on 
the underlying mass distribution of (sub)haloes, appear to 
be largely model-independent. To the precision achievable 
with our simulations, the failures previously attributed to 
the CDM paradigm appear to be largely attributable to a 
neglect of baryons in dark matter only simulations. 

Overall, we believe that the GIMIC simulation, which 
has already been demonstrated to s uccessfully repro- 
duce many aspec t s of galaxy formation JCrain et al. 201C 



duce m any aspect s 01 galaxy tormation llUram et aLIUUU 
Font et al] 1201 ll . IMcCarthv et all |2012| . IMcCarthv et al 



2QLj), 



is a much closer approximation to the observable 
Universe than a dark matter only model. With current and 
upcoming surveys such as ALFALFA, GAMA and GAIA 
pushing the observations to even fainter limits, simulations 
must not only increase in resolution, but will also have to 
take baryons into account, if they are to resolve the scales 
required to distinguish alternative dark matter models, or 
to model galaxy formation on sub-Milky- Way scales. 

While it seems somewhat unsatisfactory to modify ap- 
parently "assumption-free" and elegant methods like abun- 
dance matching and to replace dark matter simulations with 
complex, and in many ways uncertain astrophysical simula- 
tions, it was clear from the outset that a "dark matter only 
universe" itself was merely an assumption of convenience. 
On small scales, it falls short. 
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APPENDIX A: MASS CORRECTION FOR 
CENTRALS AND SATELLITES 

In Section \3. 2 1 we parametrised the average relative change 
of a subhalo's mass from the DMO to the GIMIC simulation 
in Eq. [2] 

Mqimic = a + (MpMo/Mt)" 
Mdmo l + (M D Mo/M t r 

In Fig. IA1I we show the ratio of subhalo masses for matched 
pairs of satellites, centrals, and the combination of all sub- 
haloes at four different redshifts, from z = 6 to z = 0. 
Allowing a, M t and w to vary freely, we fit Eq. [2] to the 
different populations of subhaloes, giving equal weight to 
the median ratio within each logarithmic mass bin. It can 
be seen that for central subhaloes, the lower limit does not 
evolve strongly with redshift from z = 1, although the mass 
M t at which the mass ratio reaches the intermediate value 
of (a+ l)/2 rises from 10 n - 2 M Q at z = 1 to lO n - 6 M at 
z = 0. At 2 = 6, the reduction in subhalo mass is smaller. 

We noted in Section [3. ll that at fixed subhalo mass at 
z = 0, the average difference in mass between the GIMIC 
simulation and the DMO simulation is slightly less for satel- 
lites than for centrals, which may be attributed to the fact 
that satellites also experience tidal effects, which are similar 
in both simulations. A pair of isolated subhaloes that has 
evolved to a given mass ratio before infall, and whose mass 
is subsequently reduced in equal proportion by tidal effects 
in both simulations, would evolve to the same mass ratio, 
but for a lower total mass. Fig. IA1I shows that a difference 
also exists at higher redshifts, where the decrease in mass ra- 
tio appears at higher masses for satellites than for centrals, 
although at high redshift, the number of massive satellites 
is low, and the associated statistical uncertainty is large. It 
is worth noting that the fraction of satellites increases with 
time, and that more than half of the satellites at z = were 
still (isolated) centrals up to z = 0.5, partly explaining the 
convergence of the two curves over time. At all times, the 
average reduction in mass for the total population of sub- 
haloes closely resembles that for the centrals. 
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Figure Al. Ratio of subhalo masses between the GIMIC and DMO simulations for matched pairs of subhalocs at different rcdshift. The 
error bars indicate the estimate of the median and its error for centrals (blue), satellites (red) and all subhaloes (black), while the lines 
show the corresponding best fits to Eqn. [2] with coefficients listed in Table lAll 
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logioMt 


w 
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0.40 
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0.78 


11.07 


5 


6.0 


all 


0.77 


11.18 


1.95 




centrals 


0.77 


11.20 


1.58 




satellites 









Table Al. Coefficients of Eq. [2]for fits to the median mass ratios 
of individual, matched subhaloes in the GIMIC relative to the 
DMO simulation. For each redshift, the rows show the result of 
different sets of subhaloes, as defined in Section 12.2.31 Note that 
at high redshift, scatter dominates the mass ratio for satellites, 
but the total population remains dominated by centrals. 



Table IA1I contains the coefficients of the different fits, 
and may be useful for constructing subhalo mass functions 
based on a DMO catalogue. 



APPENDIX B: CONVERGENCE 

As described in Section [2] the GIMIC and DMO simula- 
tions are each performed at two different resolutions, with 
particle masses that differ by a factor of 8. In Fig. IB1I we 
compare the cumulative mass functions in both sets of sim- 
ulations, with thin and thick lines denoting the low- and 
high- resolution results respectively. The left panel shows 
the halo mass function, while the right panel shows the sub- 
halo mass function, both at z = 0. In both panels, black 
lines denote the DMO simulations and red lines the total 
mass functions of the GIMIC simulations. While the total 
number of both haloes and subhaloes decreases by an ex- 
pected factor of ~ 5 for the lower resolution runs, it can 
be seen that the mass functions are well converged, up to 
the resolution limit. Hence, the baryonic effects on the total 
number of haloes and subhaloes as a function of mass, as 
measured in the GIMIC simulation, are largely independent 
of resolution, at least for masses of ~ 10 9 Mq and above. 

In the right panel of Fig. IBU we also plot the reduced 
cumulative mass functions of subhaloes that contain stars 
(green), or gas (blue). While the reduced mass functions for 
subhaloes with gas are well converged, the reduced mass 
functions of subhaloes with stars begin to diverge above the 
absolute resolution limit, indicating some resolution depen- 
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Figure Bl. Cumulative mass function of haloes (left panel) and subhaloes (right panel) at z = in the high-resolution (thick lines) 
and intermediate resolution (thin lines) realisations. In both panels, results from the DMO simulations are shown in black, while those 
of the GIMIC simulations arc shown in red. Also plotted in the right panel are the reduced subhalo mass functions for subhaloes with 
stars (green), and with gas (blue). It can be seen that the total number of haloes and subhaloes are converged to the resolution limit, 
both in the DMO and the GIMIC simulations. While the reduced mass functions for subhaloes with gas are also converged, the reduced 
mass functions of subhaloes with stars begin to diverge at higher masses, indicating some resolution dependence. 

dence of the star formation threshold. In addition to model 
dependence, this lack of convergence implies some uncer- 
tainty in the effect on stellar-to-total mass ratios discussed 
in Section [4] but does not change our finding that a sig- 
nificant fraction of dark subhaloes are present in a universe 
with baryons, which has a significant effect on abundance 
matching at the low-mass end. 
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